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ABSTRACT 

As India has pushed into and beneath the south margin of Asia in Cenozoic time, 
it has added a great volume of crust, which may have been either (i) emplaced 
locally beneath Tibet, (ii) distributed as regional crustal thickening of Asia, (Hi) 
converted to mantle eclogite by high-pressure metamorphism, or (iv) extruded 
eastward to increase the area of Asia. The amount of eastward extrusion or 
"escape" is especially controversial: plane-stress computer models of finite strain 
in a continuum lithosphere have shown minimal escape, while laboratory and 
theoretical plane-strain models of finite strain in a faulted lithosphere have shown 
escape as the dominant mode. We suggest that the best way to work toward a 
resolution is to study the present (or «eo)tectonics by a set of computational 
experiments, making use of the known fault network and using available data on 
fault activity, geodesy, and stress to select the best model. 

We apply a new thin-shell method which can represent faulted lithosphere 
of realistic rheology on a sphere, and provide predictions of present velocities, 
fault slip-rates, and stresses for various trial rheologies and boundary conditions. 
To minimize artificial boundaries, the models include all of Asia east of 40°E and 
span 1 00° on the globe. The primary unknowns are the friction coefficient of 
faults within Asia and the amounts of shear traction applied to Asia in the 
Himalayan and oceanic subduction zones around its margins. Data on Quaternary 
fault activity prove to be most useful in rating the models. 

The best results are obtained with very low fault friction of 0.085 (only 
10% of the friction assumed within unfaulted blocks): this major heterogeneity 
shows that unfaulted continuum models cannot be expected to give accurate 
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simulations of the orogeny. On the other hand, even with such weak faults, only a 
fraction of the internal deformation is expressed as fault slip; this means that rigid 
microplate models cannot represent the kinematics either. A universal feature of 
the better models is that eastern China and southeast Asia flow rapidly eastward 
(e g., 20-60 mm/a) with respect to Siberia. The rate of escape is very sensitive to 
the level of shear traction in the Pacific subduction zones, which is below 6 MPa. 
Because this flow occurs across a wide range of latitudes, the net eastward escape 
is greater than the rate of crustal addition in the Himalaya. The crustal budget is 
balanced by extension and thinning, primarily within the Tibetan plateau and the 
Baikal rift. The low level of deviatoric stresses in the best models suggests that 
topographic stress plays a major role in the orogeny; thus, we have to expect that 
different topography in the past may have been linked with fundamentally 
different modes of continental collision. 

INTRODUCTION 

Our zeroth-order understanding of why continental collisions (orogenies) differ 
from normal oceanic subduction is that continental crust is too buoyant to 
subduct, and accumulates near the surface in a thick layer which supports high 
topography and thus changes the patterns of tectonic stress and flow. The 
beginning of this process in the Himalaya is well dated at about 55 Ma (see 
chapters by Rowley and Burbank in this volume). Using the known plate motions 
since this time, the length of the Himalayan arc, and the crustal thickness of India, 
it is relatively straightforward to compute that about 4± 1 x 10 5 km 3 of continental 
crust have been intruded into or subducted under Asia since that time (Le Pichon, 
Fournier, and Jolivct. 1992). The most basic question one can ask about the 
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orogeny is, where has this excess crust gone? 

Fate of the excess crustal volume 

For many years, the debate about this question has emphasized two opposite 
ideas: that the excess crust was stacked locally beneath Tibet, or was extruded 
eastward to reduce the area of the Pacific basin. 

Growth of Tibet. Since Tibet is the highest and widest plateau on the Earth, and 
coincides exactly in longitude with the Himalaya, it is clearly a result of the 
orogeny. Crustal thickness here is about 70 km, or twice that of normal 
continents. The volume of crust deposited here is uncertain by about a factor of 
two, depending upon whether the pre-orogenic topography resembled the 
Altiplano of the Andes, or the low arc of Sumatra. Also, it is very controversial 
how the crustal thickness was doubled. Emile Argand envisioned a single 
underthrust of India; this view was recently restated by Powell and Conaghan 
(1973). A recent variant of this idea is that underthrusting India may have pushed 
a "pool" of soft lower crust of Asia ahead of it, "inflating" Tibet more uniformly 
(Zhao and Morgan, 1987; Bird, 1991). Others have argued that underthrusting is 
minor, and that Tibet has been horizontally shortened to about 50% of its former 
width (Dewey and Burke, 1973). The idea that Tibet has absorbed most of the 
surplus crust is supported by a number of plane-stress, thin-plate, finite-strain 
continuum computer models (e.g., England and McKenzie, 1982; Vilotte, 
Daignieres, Madariaga, and Zienkiewicz, 1982). When such models are given 
rigid boundary conditions on all sides, eastward escape is impossible; yet even if 
escape is permitted by a free eastern boundary, it is never very important in such 
continuum models (see chapter by Houseman in this volume). In the 
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homogeneous flat sheet of lithosphere around the model "Tibetan plateau", strain- 
rates typically decay as a power of distance from the highlands, precluding 
anything resembling through-going master transform faults. Certain specific 
heterogeneities (like a rigid Tarim Basin) have been added to some of these 
models (Vilotte et al., 1984; England and Houseman, 1985) and this causes local 
perturbations and broad shear zones, but does not change the basic result. 
Continental escape. The antithesis in this debate is that China was extruded to the 
east, overriding and shrinking the Pacific basin. This idea originated in the 
mapping (from satellite images) of great transform faults such as the Altun Tagh 
and Red River faults, the linkage of Red River fault slip with the Cenozoic 
opening of the South China Sea, and a proposed analogy between Asia and a sheet 
of plastic material deformed in plane-strain (Molnarand Tapponnier, 1975; 
Tapponnier and Molnar, 1976; 1 977). (According to this view, gravity prevents 
further crustal thickening after an initial transient, so the later phases of the 
orogeny have minimal strain along vertical axes.) Transform-dominated extrusion 
tectonics of this style was modeled using strain-weakening plasticines or clays 
(which form faults spontaneously) with a plane-strain constraint (Tapponnier, 
Peltzer, Le Dain, and Armijo, 1982; Peltzer and Tapponnier, 1988). 

Actually, the debate should be widened to include two other possibilities: 

Regional crustal thickening. It is notoriously difficult to determine 
paleoelevations in the absence of marine sediments, and Cenozoic marine rocks 
are found only in eastern China. Thus, one cannot rule out distributed crustal 
thickening in Siberia, Mongolia, and western China by lesser amounts, which 
might have increased the elevation of Asia by perhaps a kilometer. 
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Conversion of gabbroic lower crust to eclogite. Many believe that the typical 
composition of the lower continental crust resembles that of gabbro. If so, then 
whenever crust is thickened and heated beyond 600-800°C, its lower layers could 
convert to the rock eclogite, primarily by the growth of garnets and the dissolution 
of feldspars (Ahrens and Schubert, 1975). This transformation would convert the 
rock from a "crustal" density of about 3000 kg nv 3 to a "mantle" density like 3300 
kg nr 3 . For tectonic purposes, this volume of material would be subtracted from 
the crust and added to the mantle. It is impossible to compute a prediction of the 
size of this effect, because the transition is broad in P/T space, the reaction 
kinetics are sluggish, and the actual composition of most lower crust is unknown. 
Budget for excess crust 

In a recent attempt to quantify and account for all of the excess continental crust 
in this orogeny, Le Pichon et a/. (1992) showed that crustal thickening (both in 
Tibet and in the rest of Asia) probably accounted for 60% of the excess volume. 
This implies that about 40% of the excess volume went either into eclogite or into 
continental escape. Unfortunately, there is no known way to quantify the 
production of eclogite, except that it is probably non-negative. This only gives an 
upper bound (40%), but no lower bound, on the total amount of continental escape 
during the orogeny. Also, the assumptions behind such calculations are 
themselves controversial; by using a lower initial elevation for Asia and a greater 
volume for eroded sediment, Rowley (this volume) reaches a much lower limit 
(5%) for the amount of continental escape. The uncertainties in such an analysis 
are quite large. 
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History of the disposition of excess crust 

Can these views be reconciled by postulating that "escape" was the dominant 
mode during only one interval of the orogeny? Theoretically, one would expect 
that Tibet rose first, gradually increasing the deviatoric stress on adjacent parts of 
Asia, until the great transforms broke through and made escape the dominant 
mode in the latter part of the orogeny. However, Harrison, Copeland, Kidd, and 
Yin (1992) have reviewed the sedimentologic and metamorphic data and showed 
that Red River fault activity (the best-documented part of the escape) was 
essentially completed before the rise of Tibet. In a further complication, Molnar 
and Tapponnier (1978) documented active normal faulting in Tibet which shows 
it has actually been collapsing vertically and extending in the East-West direction 
since the Pliocene, implying that escape should once more be dominant today. 
However, it is necessary to recall the alternative possibilities of regional crustal 
thickening or eclogite formation. Hence, there is no easy solution of this kind. 
Research plan 

Considering how many fine scientists have argued over this problem for two 
decades, we cannot expect to resolve it with one new argument or computation. 
Our view is that we should begin by understanding the dynamics (i.e., the 
kinematics or flow field, and also the balance of forces) in the present, and only 
then consider how to extrapolate into the past. If at first we only try to model 
present velocities and stresses, we are freed from the need to choose initial 
conditions for the shape of Asia, the shape of India, the map of crustal thickness, 
and the map of heat-flow. Instead, we can work with published maps of present 
elevation, heat flow, and crustal thickness which are based on actual data 
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(however sparse). We know the present relative velocities of the surrounding 
plates with good precision. We can incorporate existing maps of active 
Quaternary faults, and test whether these planes are major weaknesses, or only 
passive strain-orientation indicators analogous to "slip lines" in the plastic 
deformation of metals. Especially, we can benefit from historical data on 
seismicity, geodesy, and stress directions; these will not be input data to the 
models, but instead can be compared to predicted outputs, as a test of the relative 
realism of different models in a set. 

One final motivation for beginning with neotectonic studies is that a large 
fraction of the world's population lives in Asia under the constant threat of 
devastating earthquakes. Since economic resources are not sufficient to move all 
of these people into earthquake-resistant housing, some prioritization based on 
relative risk is clearly needed. A successful neotectonic model could predict the 
long-term mean slip rates of all the major faults in Asia, which could provide 
some measure of their relative seismic hazard, until it becomes possible to replace 
these estimates with the results of local geologic and geodetic studies. 

FINITE-ELEMENT MODEL 
Assumptions 

Spherical planet. We ignore planetary rotation and approximate the lithosphere as 
a spherical shell (of variable thickness). We assume radial, laterally- 
homogeneous gravity, so that the sphere is a geoid. (As with mantle-convection 
theories, our methods can be applied to a real rotating planet by the simple trick of 
measuring relative radial position with respect to the actual geoid, not with respect 
to the planet's center.) 
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Creeping flow. Because we seek to compute velocities which are averaged over a 
longer time scale than that of the earthquake cycle, we ignore individual 
earthquakes and all other accelerations except gravity. 

Anelastic rheology. In a quasi-steady state with time-invariant boundary 
conditions, elasticity contributes a negligible fraction of the strain-rate in 
viscoelastic solutions. The "viscous" part of the strain-rate is caused by 
thermally-activated dislocation creep and/or frictional sliding on faults, both of 
which are nonlinear ("non-Newtonian"). We neglect elastic strain entirely to 
eliminate the need for unknown initial conditions, and for time-steps. 

No flexural strength. We assume no vertical shear traction on vertical planes, and 
that vertical normal stress is therefore lithostatic at all points. (Also, most parts of 
most models will be designed to be isostatic, but this is not a constraint of the 
method. If the elevation and density structure taken together would imply 
anomalous tractions at the base of the model, such tractions will be properly 
considered in the momentum balance.) 

"Thin-plate" (actually, "thin-shell"). The horizontal components of the 
momentum equation are vertically (radially) integrated through the lithosphere 
before they are applied as constraints on the solution. In this paper, we shall also 
assume that the horizontal components of velocity are nearly independent of depth 
(actually, they are proportional to r). This second assumption is not essential, and 
Bird (1989) showed how it could be relaxed in the case of flat-Earth models. 
Constant thermal properties. We assume constant thermal conductivity and heat 
productivity in all parts of crust, with (distinct) constant properties in all mantle 
lithosphere. 
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Incompressibility (consistent with neglect of elastic strain). 

Steady state, vertical heat conduction. The geotherm is assumed to depend only 
on heat flow and thermal properties. (This assumption can also be relaxed, with 
additional programming.) 

The equations derived from these assumptions are contained in Kong and 
Bird (in press), together with a report on verification tests of the resulting 
Fortran77 program (SHELLS). Essentially, we divide the lithosphere into 
spherical-triangle elements, whose corner points are called nodes. The strength of 
the lithosphere is determined by 3-D numerical integrals throughout the volume, 
taking account of different crust and mantle properties, varying geotherms, 
varying strain-rates, and migration of the brittle/ductile transition(s). The 
horizontal components of velocity at nodes are the primary unknowns that we 
solve for. Because the rheologies we use are nonlinear, it is necessary to 
approximate them with linear ones, and then iterate each computation until the 
linearized and nonlinear rheologies coincide at all integration points in the grid. 
Input Datasets 

Grid geometry. The finite element grid used in our simulation experiments is 
shown in Figure 1 . All previous simulations of Asia (whether in the computer or 
in the lab) have used plane methods. This made it necessary either to move the 
model boundaries artificially close to Tibet, or to accept geometric distortions 
near the model edges. In practice, all previous authors have used artificial straight 
boundaries. We take advantage of the spherical-shell method to use mostly 
natural plate boundaries for Asia. Along the south margin, we include the Zagros 
orogen, the Makran oceanic subduction zone, the sinistral transforms of Pakistan, 
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the Himalayan frontal thrust, the dextral transpressive boundary of Burma, and the 
Indonesian subduction zone. (In each case, the edge of the grid is a fault 
element.) On the east, the boundary follows the oceanic subduction zones of the 
Philippines, Taiwan, Japan, and the Kamchatka arc. In the northeast, we follow 
the North America/Eurasia boundary proposed by Cook, Fujita, and McMullen, 
(1986) to the Nansen Ridge. It is only at this point that we depart from actual 
plate boundaries, to avoid modeling the complex Mediterranean region and the 
eastern Atlantic. Our model is terminated at the meridian of 40°E (west of the 
Urals), beyond which the Himalayan orogeny is unlikely to have any influence. 

The faults included in the grid include all those indicated as active or 
potentially active on two tectonic maps of Asia: Institute of Geology (1979) and 
Ma (1987). We also include a single thrust fault to represent the inactive Ural 
suture, to see if the stress fields in our models would disturb it. Because it is 
unlikely that mature faults of large displacement actually terminate without 
connections, we connect together adjacent faults in places where young sediment 
may prevent the connection from being mapped. (This connectivity is quite 
important to preserve the effective degrees of freedom of models in which the 
faults are weak.) 

It is necessary to assume a dip for each fault segment. We first classify all 
faults tentatively as thrust, normal, or strike-slip based on field data, fault plane 
solutions, topography, or mechanical intuition. We then assign dips of 25° for 
thrusts, 65° for normal faults, and 90° for strike-slip faults, consistent with the 
hypothesis that they initially formed as Mohr-Coulomb shear fractures in a 
homogeneous medium of friction / = 0.85. It is important to realize that the 
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sense of dip-slip tentatively assigned is not a constraint on the solutions: a fault 
assigned a dip of 25° may actually slip in a thrust, transpressive, strike-slip, 
transtensional, or normal sense, depending on the stresses applied to it. Only 
faults assigned a dip of 90° are constrained to move in a strike-slip sense, either 
dextral or sinistral. (To avoid artificial locking at a joint between two strike-slip 
fault elements of different azimuths, their mean azimuth is used to impose this 
constraint.) Because an incorrect strike-slip constraint could eliminate necessary 
degrees of freedom from our solutions, we assign a compromise dip of 45° to 
some faults whose dip and sense of slip are uncertain (e g., the boundary fault 
between India and Burma). All faults will lock if there is not enough deviatoric 
stress in their vicinity. 

Elevation. To accurately represent the forces associated with topography, 
elevations and ocean depths are taken from the ETOP05 global dataset of 5' 
means. Then, elevations of nodes are adjusted so that the piecewise-planar 
topography of the model is best-fit to this data. Finally, the depths of all 
subduction-zone trenches are corrected by hand, since the process above tends to 
underestimate them. Accurate trench depths are important so that the lithostatic 
component of the boundary tractions on the model edges will be correct. 
Heat-flow. Unfortunately, there are large parts of Asia with no heat-flow wells, at 
least in the open literature. Our map of heat-flow (Figure 2) is based on 5° means 
kindly provided by Henry Pollack (personal communication, 1993), with some 
liner-scale detail based on maps by Duchkov et al. (1985) (Siberia) and Lysak 
( 1 992) (Baikal rift). Excluding the area around the spreading Nansen Ridge, this 
heat-flow model is conservative, ranging only from 52 to 100 mW nr 2 . Most of 
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Tibet and the tectonically active area to the east (as far as 1 10°E) is represented by 
a uniform 60-67 mW nr 2 . 

Lithosphere thickness. We assume that the base of the mantle lithosphere is an 
isothermal surface at 1300°C, and compute its depth from the heat-flow data, 
assuming thermal conductivities of 3.5 and 6.0 W nr 1 “K' 1 in the crust and 
mantle, respectively, and radioactive heat production of 4.6x1 O' 7 W nr 3 in the 
crust. The thickness of the mantle lithosphere is typically 85-105 km, but values 
up to 160 km are inferred in the Siberian craton, and smaller values are computed 
for East Tibet, the Sea of Japan, and in the Baikal rift (Zorin, Kozhevnikov, 
Novoselova, and Turutanov, 1989). 

Crustal thickness. After the computations above, we adjust the crustal thickness 
to achieve local isostasy everywhere with respect to a typical midocean rise. This 
computation assumes crust and mantle densities of 2875 and 3330 kg nr 3 , 
respectively (at 0°K; the corresponding volumetric thermal expansion coefficients 
are 2.4 and 3.1x1 0 -5 K 1 ). The resulting thicknesses range from 5 km (enforced as 
a minimum) in the ocean basins to 75 km in Tibet, with a modal value of 40 km in 
Asia. Because of the constraint that crustal thickness is no less than 5 km, a very 
thin strip of model along each oceanic trench can not be isostatically balanced; it 
is assumed that the necessary tensional anomaly in the vertical stress at the base of 
the lithosphere is provided by the dynamics of subduction. 

Rheology. All simulations share a common pair of dislocation-creep rheologies 
(one valid in the crust, and the other in the mantle lithosphere). The equation for 
creep strength involves principal stresses a, (with tension positive), principal 
strain-rales <?,, absolute temperature T , depth z, and four constants ( A , B , C, n): 
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_ (cr 3 - cr, ) < 2 A(e 3 - e , ^2^-e x e 2 - e 2 e 3 - e 3 e, f 1 *exp 

The crust/mantle values of each constant are: A = 2.3xl0 9 /5.4xl0 4 Pa s ,/n ; B = 
4,000/18,314 °K; C = 0/0.017 °K nr 1 ; and n = 3. These values were based on 
previous experience in modeling California (Bird and Kong, 1994) and Alaska 
(Bird, in press) with similar programs. 

At low temperatures, the strength of the lithosphere is limited by frictional 
sliding. Typical laboratory tests show friction coefficients f = 0.60 to 0.85 for a 
wide variety of rocks. We assume that the continuum blocks between faults have 
this friction, and that they are pre-fractured on a variety of planes (we ignore 
cohesion). However, most models assume a lower friction coefficient in fault 
elements (/ f > 0.04). All frictional strengths are computed assuming locally 
hydrostatic pore pressure. If major faults are actually filled with anomalous, near- 
lithostatic pore pressure, then their true friction coefficients must be higher than 
the apparent coefficients determined in this study. 

The many active subduction zones which ring the model on the southern 
and eastern sides (including the Himalayan continental subduction zone) are 
probably different from other thrust faults because of the rapid subduction of 
sediments and the dynamic maintenance of high internal pore pressures. In 
previous studies. Bird (1978a) found a maximum mean traction of only 20 MPa in 
the Himalayan zone, and Bird (1978b) found average shear tractions of only 9-32 
MPa in the Tonga and Mariana oceanic subduction zones. To represent this, we 
have added special programming which monitors the mean shear traction in these 
boundary subduction zones and prevents these tractions from exceeding a preset 
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limit. 



15 


Boundary conditions. Around most of the perimeter of the model the edge of the 
grid is ringed with fault elements. Therefore the outermost nodes belong to the 
adjacent plates, and it is these nodes which require boundary conditions. We 
impose the present velocities of the Arabia, India, Australia, Pacific, and North 
America plates (with respect to Eurasia) from the NUVEL-1 model of DeMets, 
Gordon, Argus, and Stein (1990). The present velocities of the Philippine plate 
are from Seno et al., 1987. The short northern boundary is held fixed in the 
reference frame of western Eurasia (Europe). Along the artificial western 
boundary, we permit free slip in the N-S direction, but prevent any E-W velocity 
component. (This particular boundary location was chosen because it lies along a 
local mirror plane of symmetry in the orocline of the Zagros orogen. Local 
symmetry then dictates this choice of boundary conditions.) 

Scoring 

In order to choose the best model, we compare our modeling results with geologic 
and geophysical data, such as relative velocities between pairs of points measured 
by Very Long Baseline Interferometry (VLBI) or Global Positioning System 
(GPS) geodesy, long-term-average fault slip rates from geological or 
morphological research, and stress directions. 

Maximum horizontal stress directions were compiled by Zoback (1992) 
for the whole Earth. There are 1,624 measurement points in our model region. In 
North China, there are many northeast-trending active normal faults, and the 
maximum horizontal stress (g, (I ) directions associated with a known normal- 
faulting stress regime are mostly NE or ENE, as we would expect. However, 
some o /H data of "unknown stress regime" point N W. Therefore, we reject o,,, 
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directions of "unknown stress regime" as less reliable and use only the 1,033 ct ai 
directions with known stress regime. From a variogram which we computed 
using this data, we find that the stress directions contain random errors with 
standard deviation of 30°. With 30° error, this dataset is not a very reliable scoring 
reference. 

Heki and Koyama (1993) reprocessed the VLBI data of Ryan et al. (1993) 
and found that Shanghai in east China is moving southeastward at 10 mm/a 
relative the interior of the Eurasia. But Robaudo and Harrison (1993) used 
combined VLBI and Satellite Laser Ranging (SLR) data and got the opposite 
result for the same station, that Shanghai is moving northwestward at 10±2 mm/a 
with respect to Eurasia. Unfortunately, Shanghai is the only station in Eurasia 
that is east of the Urals. At the present time, we cannot tell which of above two 
opposing conclusions is correct. So, we will not this datum in scoring our results. 
However, our model results will be able to tell which is most reasonable in terms 
of fitting with other data sets. 

Only a very few major faults in Asia have known slip rates. We use the 
slip rates of 1 3 major faults (Table 1 ) to score our model results in a second way. 
One common problem with published geologic slip rates is that it was difficult to 
determine the age of the offset feature. Some ages are merely assumed (e g., 
Peltzer, Tapponier, Gaudemer, and Meyer, 1988; Peltzer, Tapponnier, and 
Armijo, 1989). Another consideration in using these data is that the faults in our 
model represent fault zones, not individual faults; but the published geologic slip 
rates are often for individual faults. Other published slip rates (Ma, 1987) were 
determined from sums of seismic moments, and it is questionable whether these 
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short-term rates can represent the long-term average activity of the faults. It 
almost seems that there is no reliable data we can use to test our modeling results. 

However, for most faults in Asia, we know their motion senses, that is, 
relative movement directions. Although the slip rates of some major faults are 
controversial, there are few questions about their motion senses. We compile a 
dataset of fault motion senses for 19 major faults in Asia (Table 2), and use this to 
compute a third kind of score. 

Considering the problems of above data sets, we will give greater weight 
to fault motion sense them to the other data sets in choosing the best models. 
Experiments 

The first set of models is to test if the friction coefficient on faults in Asia is as 
small as that in southern California (Bird and Kong, 1994) and Alaska (Bird, in 
press). In this set of models, there is no shear traction limit on the subduction 
zones, which are treated like any other fault element. The internal friction 
coefficient of the blocks between faults is kept at 0.85. We calculate models with 
fault friction coefficients (/ f ) of 0.085, 0.17, 0.51 and 0.85. Table 3 gives the 
scores of these 4 models. The errors of fault slip rate and stress direction hardly 

vary. The percentage of faults with the correct movement sense declines with 
each increase of fault friction. At / f larger than 0.51, almost all faults are inactive 
(Figure 3). When / f is small (0.17 to 0.085), many faults are active (Figure 4) 
and about 60% of the predicted fault slip directions match the sense data. At j\ 
equals 0.085, the velocity of Shanghai (relative to Eurasia) is close to the geodetic 
result of Robaudo and Harrison (1993). The results of this set of models show 
that the major faults of Asia have low friction coefficient, as in southern 
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California and Alaska. However, the best model still has important defects. 

Many major faults are still not moving or move in the wrong direction. For 
example, the Baikal rift and Shanxi graben are supposed to be normal faults but in 
Model J9401 both of them are thrusting (Figure 4). The Urals, which should be 
inactive, are thrusting slowly. All of these errors probably result from too much 
E-W compression, caused by too much shear traction on the Pacific subduction 
margin. The Altun Tagh fault is not as active as it should be. And, slip rates 
along the plate boundaries with the India and Australia plates are probably too 
small. Obviously, even the best model in this group is not satisfactory. 

Bird (1978a) concluded that the mean shear traction along the Main 
Boundary Thrust between the India plate and Eurasia plates should not be larger 
than 20 MPa. The slow slip of this fault in model J9401 (6-7 mm/a) also implies 
that we are overestimating the strength of this fault. For our next set of 
experiments, we test the effects of varying the traction on the interplate surfaces 
of subduction zones. We will treat the boundary thrust faults as a special kind of 
fault on which shear traction will be limited to a preset value. This value will be 
an input parameter that can be adjusted to fit the scoring datasets. We vary this 
parameter ( T) from 30 down to 5 MPa. When it is more than 20 MPa, the Main 
Boundary thrust hardly moves. From seismic moments, one can infer about 36 
mm/a slip rate on this fault (Ma, 1987). Jackson and Bilham (1993) used relative 
uplift rates of 6±2 mm/a in the Himalaya to infer a convergence rate of about 10 
mm/a for the ramp portion of the boundary thrust. Although there is controversy 
about the slip rate of the Himalayan boundary, we are sure that the Main 
Boundary thrust is quite active. If we take the numbers above as the range of the 



19 


slip rate of the Main Boundary thrust, the slip rate should be 10 to 36 mm/a. 

When T equals 15 MPa, the slip rates along the Main Boundary fault fall within 
this range. When this parameter is less than 10 MPa, the Tibetan plateau 
completely collapses along its southern edge (Figure 5). A shear traction of 15 
MPa along the Main Boundary Thrust is probably a reasonable estimate, and we 
will use this value for our following calculations. 

However, the fault motions inside the Asian continent is still not right. 

The major faults on the eastern side of the Tibetan plateau are not moving. The 
Altun Tagh is too slower and Baikal is almost a purely strike-slip fault. To 
improve our modeling results our next step is to test possible differences in 
traction between the Main Boundary Thrust and the oceanic subduction zones. 

We will use two parameters to adjust the traction on subduction zone surfaces. 
One is for the traction on the Main Boundary Thrust; we call it the Himalayan 
subduction traction (T H ). The other is for all other subduction zone and is called 
the oceanic subduction traction (7^). We assign different values for both T n and 
T 0 ranging from 1 to 20 MPa. No matter how T 0 changes, as mentioned above, 
the slip rates of the Main Boundary thrust are always within the range of 4 to 36 
mm/a when T H is between 13 and 17 MPa, with an ideal value of 15 MPa. (And, 
as we will see later, with 7j, equal to 15 MPa we can get better scoring results.) 
When T 0 equals 2 to 4 MPa, we get the best scoring results for both slip rate and 

fault motion sense (Table 4). 

These scores of slip rate and fault motion data favor low shear traction on 
the oceanic subduction zones. But the stress direction errors are large (40°) in 
each model. We do not consider this a serious problem. Even if our modeling 
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results were only a few degrees different from the true stress direction, the stress 
dataset has 30° standard deviation, so it must be expected that the mismatch of 
model to data would be more than 30°. Stress is apparently not well enough 
known to discriminate between models. 

Our best model has parameters: Himalayan traction 7^ = 15 MPa, oceanic 
traction T 0 =2 MPa, fault friction / f = 0.065, and block friction = 0.85. This 
model also predicts a velocity at Shanghai, China close to that of Ryan et al. 
(1993) and Heki and Koyama (1993). 

Figure 6 shows our predicted velocity field in Asia. Material is moving 
away from the Tibetan plateau toward the eastern boundary of Asia. This result 
shows that material is currently "extruding" from the Himalyan orogen. This 
eastward motion is mainly to the SE of a line from the Pamirs to Baikal. On the 
NW side of this line the crust is very stable. This extrusion is distributed over 
almost the whole of central and east Asia, and is not the result of rigid block 
motion. A major extrusion boundary lies along the Tian Shan, instead of at the 
Altun Tagh fault. 

Unfortunately, we do not know when this kind of extrusion started, 
whether at the collision between India and Asia, or after the initial rise of Tibet, or 
only after the postulated delamination of Tibet in the Pliocene. Our modeling 
cannot address such historical questions until palinspastic maps of 
paleotopography and fault networks have been drafted. 

Figure 7 shows our predictions of long-term average fault slip rates. Our 
sense predictions are correct for most major faults in Asia; the Baikal rift is 
opening, and the Altun Tagh, Red River, and South and North Tian Shan faults 
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are moving in the correct direction. Our prediction of the Altun Tagh fault rate 
(2-8 mm/a) is much smaller than that of Peltzer et al. (1989), but is almost the 
same as that determined from offset paleo-stream-channels with l4 C ages (Altun 
Tagh Active Fault Research Group, 1 992). Based on a lot of field work, they 
found the dextral slip rate of the Altun Tagh fault is 4.3 to 5.2 mm/a, and the 
vertical slip rate is a minimal 0.3 to 0.8 mm/a. If it is true that Altun Tagh fault 
does not move very fast, our prediction would be that the boundary of the 
extruding material lies along the Tian Shan. 

In this model the normal faults representing the Baikal rift slips at up to 35 
mm/a, implying up to 15 mm/a of spreading. (Note: To convert from dip-slip 
rates to discontinuities in horizontal velocity, reduce thrusting rates by 9% and 
normal-faulting rates by 58%. ) This is almost equal to the relative velocity of the 
surrounding blocks, which is about 1 8 mm/a. Thus, the model suggests that the 
rifting is primarily "passive" (driven by microplate motions) rather than "active" 
(driven by local density anomalies). The main role of the thin lithosphere and 
high heat-flow in the Baikal rift is to cause the spreading to concentrate there, 
whereas to the SW it is very diffuse (Figure 8). 

One big problem with our results is that the Xianshui He fault and the 
Longman Shan thrust are not moving. This might result from the low resolution 
of our model, especially the very smooth heat flow field which we used (Figure 
2). In the region of these faults, almost no heat flow data is available. So heat 
flow is uniform across these faults in our model, and the rheology (which is 
dependent on the geotherm) is uniform. Perhaps some special local weakness is 
required to facilitate the slumping of crust down the eastern margin of Tibet. 
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Because of this failure in our prediction, at this stage our model is not yet ready to 
estimate the seismic hazard distribution within Asia. 

Figure 8 shows the predicted principle strain-rates within the continuum 
blocks (ignoring slip on fault elements). This can be interpreted as the 
deformation due to all faults too small to appear in our grid. The deformation is 
almost all tensional or transtensional and is concentrated in the Tibetan plateau 
and its the neighboring Baikal rift and Ordos plateau. (There is thrusting in this 
model, but it occurs in fault elements where we have lowered the friction, not as 
distributed deformation in the stronger blocks.) The Tarim basin and southeast 
China are stable. 

The vertically integrated stress anomalies are shown in Figure 9. The 
quantity displayed is the vertical integral through the lithosphere of the stress 
anomaly (relative to the center of a midocean rise). The vertical component of the 
integrated stress anomaly is dominant in the Tibetan and Ordos plateaus. In the 
Tian Shan and other areas near to Tibet the maximum compressional components 
are in directions radiating from the plateau. The distribution of the stresses means 
that the weight of the plateau is very important in driving the neotectonics of Asia. 

CONCLUSIONS 

Unfortunately, even the best of our models still fail to reproduce the right 
direction of slip on a few (11% oF) active faults. In other cases, the sense is 
correct but the rate errors are large. Thus, it would be premature to use these 
models for seismic hazard estimation in Asia. To achieve major improvements 
will require more data, of which the most important would be additional heat-flow 
values for input to the mode's, and additional geologic slip-rates or geodetic 
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velocities for testing of the output. 

Even at this early stage, it is apparent that the faults of Asia share the 
anomalous weakness that was previously found in California and in Alaska. 
Previous continuum finite-strain simulations of this orogeny with finite elements 
have been based on an implicit assumption that faults are passive strain markers, 
but not important heterogeneities; this seems to be wrong. Since rapid fault slip 
implies local concentrations of straining in the ductile part of the lithosphere, 
there must be a compensating weakness in the frictional part of the fault, with 
respect to adjacent blocks. The contrasts that we find are too great to ascribe to 
variations in lithology. More likely, they are due to concentrations of pore 
pressure along faults, where previously subducted water may be returning to the 
surface. 

The thrust faults of the marginal subduction zones are even more profound 
heterogeneities: in our model, the best strengths are 15 MPa for the Himalayan 
thrust, and about 2 MPa for the oceanic subduction zones. Such low strength 
almost certainly requires near-lithostatic pore pressure in the subducting 
sediments. The pore pressure may actually equal the minimum compressive 
stress (a 3 ) in these zones, and water may escape by hydrofracturing. Or, pressure 
may be self-regulating, by a mechanism in which permeability increases 
exponentially as pore pressure approaches a 3 . In the latter case, it would be 
natural for the equilibrium pore pressure (and thus the fault strength) to depend on 
the permeability of the subducted sediment. The strength difference we infer 
between the continental and oceanic subduction zones might be explained by the 
porosity difference between sandstones and abyssal cherts, for example. 
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An approach to neotectonics that geologists often use is the geometric 
construction of microplate kinematic models, in which all the deformation is 
expressed as fault slip ( e.g ., Bird and Rosenstock, 1984). If our results are 
correct, then this approach to Asia is also faulty. We find that only a fraction of 
the extrusion in our best models can be accounted for by slip on fault elements; 
much more results from distributed transtensional straining of the blocks between. 
To date, the only modeling method that seems capable of including weak faults, 
distributed strain, and also finite strain is the laboratory deformation of scaled 
clay-cakes and other strain-weakening materials. 

Our contribution to the debate about the historical importance of eastward 
extrusion in the Himalayan orogeny is a clear result that it is very rapid today. 

We cannot state a precise rate because in our models this is sensitive to changes of 
as little as 1 MPa in the shear traction of oceanic subduction zones along the 
Pacific margin. However, it is clear that all the best models show rapid extrusion, 
and that the rate of crustal escape to the east currently exceeds the rate of resupply 
along the Himalaya. The budget is balanced by distributed crustal thinning and 
extension, which quantitatively exceeds the internal shortening by thrusting in the 
Pamires and Tian Shan. 

Because of the weak subduction boundaries and weak internal faults, the 
stresses are quite low in our preferred models of Asia (Figure 9). There are only a 
few places where the vertically integrated difference between principal stresses 
exceeds lxlO 13 N/m (i.e., 100 MPa x 100 km). Therefore, the regional stress 
field is naturally dominated by the compression radiating to Tibet, whose gravity- 
density moment is of comparable size. The present topography largely determines 
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the stress field, which determines the strain patterns and thus the evolution of 
topography. The conclusion that we draw from this is that the solution was 
probably very different (and extrusion much less important) in previous times 
when the plateau was not so high. 
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Table 1 . Slip rates of some major faults in Asia 


Fault 

Strike-slip rate (mm/a; 
dextral positive) 

Relative vertical iate 
(mm/a, thrust positive) 

Reference 


Minimum Maximum 

Minimum 

Maximum 


Altay 

4.6 

10 

0.2 

1.3 

Ma, 1987 

Altun Tagh 

-30 

- 4.3 



Peltzer et al 1 989; Altun 
Tagh Active Fault 
Research Group, 1 992 

E. West Kunlun 



1.5 

7.5 

Ma, 1987; Avouac and 
Peltzer, 1993 

Haiyan 

-10 

-5 

0.5 

0.8 

Ma, 1987 

Qilian Shan 

- 7.7 

-3 

0.4 

1.9 

Ma, 1987; Peltzer et a/., 
1988; Zhang et al 1991 

Red River 

2 

5 



Allen et al 1984 

Shanxi 

0.5 

2 


- 0.5 

Ma, 1987 

Tanlu 

1 

2 



Ma, 1987 

W. West Kunlun 

10 

25 

0.5 


Ma, 1987 

Weihe 

-6 

-2 



Ma, 1987 

Xianshui He 

-20 

-10 



Allen etai, 1991; Ma 1987 

Xiaojiang 

-7 

- 4.6 



Chen and Li, 1988; Ma 
1987 

Yinchuan 

3 

7 

- 1.0 

- 0.5 

Ma, 1987 



Table 2. The motion senses of some major faults in Asia 
(from Ma, 1987; Institute of Geology, 1979) 


Fault 

Strike-slip (D = 
dextral, S= 
sinistral) 

Dip-slip (T = 
thrust, N = 
normal) 

Afghanistan-T urkey 

D 

T 

Altay 

D 


Altun Tagh 

S 


Baikal 

S 

N 

Heto 

s 

N 

Longman Shan 

D 

T 

North Tian Shan 


T 

Qilian Shan 

s 

T 

Red River 

D 

N 

Sagaing 

D 


Shanxi 

D 

N 

South Tian Shan 

S 

T 

Tanlu 

D 

T 

Weihe 

S 

N 

West Kunlun 

D 

T 

Xianshui He 

S 


Xiaojiang 

s 

T 

Yinchuan 

D 

N 

Ural 

locked 

locked 
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Table 3. Scores of models with different fault friction 


Mode! 

Fault friction 
coefficient ff 

RMS prediction error 
of fault slip rate 
(mm/a) 

Percentage of fault slip 
senses matching data 

Average prediction 
error of stress 
direction 

J9401 

0.085 

4.31 

63.9% 

39.0° 

J9402 

0.170 

4.15 

59.6% 

38.8° 

J9403 

0.510 

3.96 

48.7% 

38.7° 

J9404 

0.850 

3.99 

42.4% 

38.6° 
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Table 4. Effects of varying shear traction of oceanic subduction 
(T 0 ) with T h = 15 MPa 


Model 

Oceanic 
shear traction 
(T 0 ), MPa 

RMS Prediction 
error of fault slip 
rate (mm/a) 

Percentage of fault 
motion sense which 
match data 

Average 
prediction error 
of stress 
direction 

19401 

2 

3.04 

88.6% 

39.7° 

19402 

4 

3.32 

88.6% 

39.3° 

19403 

6 

3.48 

84.8% 

39.7° 

19404 

8 

3.60 

77.8% 

40.3° 

19405 

10 

3.72 

75.3% 

40.0° 
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FIGURE CAPTIONS 

Figure 1. Finite element grid used for all neotectonic models. (As in all figures, a 
stereographic projection of data on the sphere of the Earth.) Solid lines are fault 
elements, with different tic marks to indicate dip: none = 90°, straight = 65°, box 
= 45°, triangle = 25°. Dashed lines are boundaries of continuum elements; these 
have no mechanical significance, but are shown to indicate resolution. The model 
domain in the entire Asian lithosphere, or all of the Eurasia plate east of 40°E. (A 
small corner of the Indian plate in the Assam syntaxis is also included because 
seismicity suggests that this area is deforming.) 

Figure 2. Assumed distribution of heat-flow, used to compute lithosphere 
thickness and strength. The unshaded contours above 120 mW nr 2 are only used 
along the Nansen Ridge. Most of the tectonically active region is modeled with 
heat-flow of 60-70 mW m' 2 . References in text. 

Figure 3. Fault slip rates from an unsuccessful model (J9403), in which the fault 
friction is too high ( / f = 0.5 1). Locked faults are dashed (except that all 
boundary faults are solid even if locked). Active faults are labeled with the 
relative slip rate in the fault plane, in mm/a. (To find the discontinuity in 
horizontal velocity, reduce thrusting rates by 9%, and normal-faulting rates by 
58%. ) Very small arrows in hanging walls show the directions of slip. In this 
model the Himalayan thrusts would be inactive, as well as the Indonesian 
subduction zone and the Altun Tagh fault. Relative plate motion would be 
expressed entirely as distributed permanent straining. 

Figure 4. Fault slip rates of a more successful model (J9401). Conventions as in 
Figure 3. Relative to the model in Figure 3, the fault friction f\ has been reduced 
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to 0.085. This causes all important faults to become active, although many move 
too slowly. Excessive E-W compression from the Pacific subduction zone margin 
causes thrusting in the Urals and in the Baikal rift and in Shanxi, all of which are 
incorrect. 

Figure 5. Close-up of predicted velocities in the Himalayan/Tibetan region, from 
a model (F9402) which fails because the assumed shear traction on the Himalayan 
thrusts is too low (10 MPa). The eastern front of the Himalaya collapses in a 
great landslide traveling over 250 mm/a south with respect to Eurasia (or 300 
mm/a with respect to India). 

Figure 6. Predicted velocities of the best model found to date (K9401). This 
model has friction of / f = 0.065 in faults but 0.85 in intervening blocks, with 
shear tractions of 1 5 MPa in the Himalaya and 2 MPa in the Pacific margin. Note 
the eastward extrusion of much of China at about 40 mm/a. (This rate is very 
sensitive to the level of oceanic subduction traction, but the pattern is not.) 

Figure 7. Fault slip rates from the model of Figure 6. Conventions as in Figure 3. 
Note 19-44 mm/a of horizontal convergence at the Himalayan front, 2-8 mm/a 
dextral slip on the Altun Tagh fault, up to 13 mm/a convergence in the Tian Shan 
and 1 8-27 mm/a in the Pamirs, up to 15 mm/a of divergence at the Baikal rift, and 
1-3 mm/a of transtensional faulting in Shanxi. Because of distributed 
deformation, these slip rates are not compatible with a rigid microplate 
description of the orogeny. 

Figure 8. Close-up detail of distributed straining (excluding faults) from the 
model of Figures 6-7. The strain-rate tensor is expressed in terms of equivalent 
faults. The many rectangular symbols indicate grabens. The "X" overprinted on 
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some indicated a component of strike-slip faulting, with the extension axis in the 
obtuse angle. While this model does include thrust faulting, it usually occurs on 
weak fault elements (see Figure 7) and not in the stronger intervening blocks 
(shown here). The only exception is some distributed thrusting (dumbell symbol) 
in the Hindu Kush. 

Figure 9. Vertically-integrated (through the lithosphere) stress anomaly (with 
respect to a midocean rise) tensors from the model of Figures 6-8. Circles 
indicate a compressive component along the vertical axis, and are proportional to 
the gravity/density moment (almost proportional to topography). In weak areas 
like Tibet, the horizontal components tend to become equal, causing horizontal 
compression to radiate through the lowlands around Tibet. 
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